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ABSTRACT 

The Near Infrared Background (NIRB) is one of a few methods that can be 
used to observe the redshifted light from early stars at a redshift of six and above, 
and thus it is imperative to understand the significance of any detection or non- 
detection of the NIRB. Fluctuations of the NIRB can provide information on the 
first structures, such as halos and their surrounding ionized regions in the Inter 
Galactic Medium (IGM). We combine, for the first time, iV-body simulations, 
radiative transfer code, and analytic calculations of luminosity of early structures 
to predict the angular power spectrum (C7) of fluctuations in the NIRB. We study, 
in detail, the effects of various assumptions about the stellar mass, the initial 
mass spectrum of stars, metallicity, the star formation efficiency (/*), the escape 
fraction of ionizing photons (/ esc ), and the star formation timescale (t$F), on the 
amplitude as well as the shape of C\. The power spectrum of NIRB fluctuations 
is maximized when /* is the largest (as Ci oc / 2 ) and / esc is the smallest (as 
more nebular emission is produced within halos). A significant uncertainty in 
the predicted amplitude of C\ exists due to our lack of knowledge of tsF of these 
early populations of galaxies, which is equivalent to our lack of knowledge of 
the mass-to-light ratio of these sources. We do not see a turnover in the NIRB 
angular power spectrum of the halo contribution, which was claimed to exist in 
the literature, and explain this as the effect of high levels of non-linear bias that 
was ignored in the previous calculations. This is partly due to our choice of the 
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minimum mass of halos contributing to NIRB (~ 2 x 10 9 M Q ), and a smaller 
minimum mass, which has a smaller non-linear bias, may still exhibit a turn 
over. Therefore, our results suggest that both the amplitude and shape of the 
NIRB power spectrum provide important information regarding the nature of 
sources contributing to the cosmic reionization. The angular power spectrum of 
the IGM, in most cases, is much smaller than the halo angular power spectrum, 
except when / esc is close to unity, tsF is longer, or the minimum redshift at which 
the star formation is occurring is high. In addition, low levels of the observed 
mean background intensity tend to rule out high values of /* > 0.2. 

Subject headings: cosmology: theory - - diffuse radiation - - galaxies: high- 
redshift — infrared: galaxies 



1. INTRODUCTION 

We have few probes of the early universe and the first few generations of stars. We 
know that stars had to form early in order to pollute the universe with metals and reion- 
ize the universe. There is evidence that the universe was reionized at around z ~ 11, 



such as from the Wilkinson Microwave Anisotropy Probe (WMAP) satellite (IKogut et al. 



20031 : ISpergel et al.ll2003l . 120071 IPage et al.ll2007l : iDunklev et al.ll2009l : Komatsu et al 



2009|). 



Stars are efficient producers of ionizing photons, so are likely candidates for the bulk of 
reionization. These ultraviolet photons at redshifts 6 < z < 30 would be redshifted 
into the near-infrared bands. Therefore, it makes sense to search for this remnant light 
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2006). Observations of the Near Infrared 



Background (NIRB) ma y indicate that there is an excess mean background above that normal 



galaxies can account for (jDwek fc Arend 



1998 



Goriian. Wright fc Chary 



2000 
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2000; 
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2005: 



Cambresy et al.ll2001tlTotani et al. 



2001 
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Kashlinskv 20051 ) . In addition, there also appears to be a peak in the N IRB spectrum 



at 1-2 /Jin, which could represent a Lyman-cutoff signature (IBock et al.ll2006l ). However, the 
interpretation of the current observational data, in particular a ccuracy of the subtract ion 
of Zodiacal light and foreground galaxies, is highly controversial ( Thompson et al. 2007aj |bh . 
Nevertheless, any detection or non-detection of this excess light could tell us properties of 
early stars. 



In addition to the mean intensity, fluctuations in the NIRB can provide an additional 
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source of information about the 


irst generations of stars (Kashlinskv & Odenwald 2000 
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2007a 


b). Fluctuations are in general easier to study than the mean intensity because 



an accurate determination of the zero point is not needed; thus, they are less sensitive 
to the imperfect subtraction of Zodiacal light. However, as the contribution to fluctua- 
tions from low redshift populations, i.e., z < 6, can confuse the signal from higher red- 
shift populations, the level of c ontamination from low redshift populations must be esti 



2007d ; iThompson et al 



with AKARI (Matsuhara et al 



iment) (IBock et al 



2006 



mated and subtracted carefully (ISullivan et al 

2007bi 



2007 



2008 



Chary. Coorav fc Sullivan! 120081 ). Upcoming measurements 



Coorav et al.l l2007t Kashlinskv et al. 



) and C IBER (the Cosmic Infrared Background Exper- 



Cooray et al.l 120091 ) may be able to put a more solid constraint on 



what fraction of the NIRB is from high redshift stars and galaxies. 

In the previous paper we have presented detailed theoretical calculations of t he spectrum 



and m etallicity /initial-mass-spectrum dependence of the mean intensity of NIRB ( Fernandez &; Komatsu 
( 120061 ). hereafter FK06). In this paper we present calculations of the power spectrum and 
metallicity/initial- mass-spectrum dependence of the NIRB fluctuations, as well as depen- 
dence on the star formation efficie ncy and the escape fraction of ionizing pho tons. While 
the previous work in the literature (jCooray et al.ll2004l ; Kashlinskv et al.l 120041 ) relied solely 
on simplified analytical arguments, we u se, for the first time, la r ge-scal e cosmological sim- 
ulations of cosmic reionization given in llliev et al.l (120061 120071 . l2008bl ). coupled with the 
analytical calculations given in FK06, to predict the power spectrum of NIRB fluctuations. 
In this way we are able to capture the contribution from ionized bubbles surrounding the 
halos, which has been ignored completely in the previous work. 



In §[2] we outline the simulations ( llliev et al.ll2006l . 120071 l2008bl ) and in §[3]we explain the 
analytic formulas we use to obtain the luminosity of the halos and the surrounding IGM. In 
§H]we present our calculation of the luminosity-density power spectrum, Pl(/c). Predictions 
for Pi{k) and the angular power spectrum of NIRB fluctuations, Ci, are presented in § |5j 
Various parameters' effects on the results are discussed in § [6j We compare our results to the 
previous literature in § [7] and to observations in § |HJ We take a look at the constraints from 
the mean NIRB in § [HI and compute the fractional anisotropy, i.e., the ratio of the power 
spectrum and the mean intensity squared, in § [TOJ We conclude in § [TTJ 
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SIMULATION 



We use simulations from llliev et al.l ( 12006k 120071 l2008bl ). which are the first truly large 
scale simulations to include radiative transfer, and are therefore ideal for predicting the 
distribution of luminosities from high redshift stellar populations. Simulations provide the 
advantage of being able to simultaneously model the distribution of halos and the density of 
the IGM, as well as the ionization front that propagates through the IGM. We combine this 
iV-body code with radiative transfer and analytic formulas for luminosity to simulate their 
luminosity-density power spectrum. 

The particula r simulation that we use in this paper is the run "f250C" in Table I of 
Iliev et al.l (l2008bl). which was run with the cosmological parameters given by the WMAP 
3-year results JSpergel et alJ l2007h . (fl m , Q A , Q h , h, a 8 , «=(0.24, 0.76, 0.042, 0.73, 0.74, 
0.95). Aside from the cosmological parameters, the only free parameter in the reionization 
simulation of this kind is the production rate of ionizing photons escaping into the IGM per 
halo. We shall come back to this parameter, called / 7 /tsF> in § 13.21 



These simulations combine a high resolution i V-body code (PMFAS T, see lMerz. Pen fc Trac 



20051 ) with a radiative transfer code (C 2 -Ray, see iMellema et al.ll2006l ). which is a conserva- 
tive, causal ray-tracing radiative transfer code. The C 2 -Ray code traces the ionization front 
by tracking photons using photon conservation. The code allows for large time steps and 
coarse grids without loss of accuracy. 

The box size of the simulation is 100 Mpc, which is large enough to sample the 
history, geometry, and statistical properties of reionization. The number of particles is 1624 3 , 
and the density field was sampled on a lattice of 3248 3 cells. The density field was then binned 
to 203 3 cells for the radiative transfer calculations. We use the 203 3 cells when we compute 
the radiation from the IGM in § 14.21 The minimum mass of the halos is 2.2 x 10 9 M Q , which 
represents dwarf galaxies. These halos have virial temperatures of 1.2 x 10 4 K, 1.8 x 10 4 K, 
and 2.6 x 10 4 K at z = 6, 10, and 15 respectively. For these halos the dominant cooling 
process is hydrogen atomic cooling. It is important to sample these dwarf galaxies, as they 
are far more numerous than larger galaxies and may provide most of the photons needed for 
reionization. 

Even though this simulation is a very powerful tool, it is important to consider its 
limitations. Halos slightly below the resolution of t his simulati o n (10 8 to 10 9 M ) may 



also be an important source for ionizing radiation, llliev et al.l (120071 ) also did a smaller 
box-size simulation [(35 h~ l Mpc) 3 ] that resolves halos down to 10 8 M , which includes 
halos that form stars as a result of atomic cooling. These smaller halos allow the ionization 
fraction to reach 50% at an earlier epoch than the simulations that only resolved down to 
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2.2 x 10 9 Mq. However, llliev et al.l ( 120071 ) found that the redshift in which reionization was 
completed remained about t he same for the 100 h" 1 Mpc and 35 h~ l Mpc simulations, due 
to the "self-regulation' : 



sec 



Iliev et al.l 120071 . for details) 



The results discussed in this paper are based on the larger box size with halos resolved 
down to 2.2 x 10 9 M & . It is possible that the smaller halos would affect the fluctuations in the 
NIRB from both the halos and the IGM. Future simulations will allow both a larger box size 
along with a smaller minimum mass. These future simulations will be able to provide more 
robust predictions for the fluctuations in the NIRB if these smaller halos contribute to the 
NIRB. Simulations that resolve halos smaller than 10 s M Q may not be needed, however: while 
these minihalos were likely the sites of the truly first generation of stars, they may not be a 
significant source of ionizing photons to reionize the universe, as UV photons in the Lyman 



(Haiman, Rees & Loeb 


1997 




Haiman. Abel & Rees 


2000 




Mac 


lacek. Brvan & Abel 


2003; 


Yoshida et al. 


2003 




Johnson, Grief & Bromm 


2008; 


Ahn et al. 


2009). However, there is 



on-going discussion as t o what the radiation feedback actually does for the formation o f 
second generation stars (IWise et al.l 120071 : lAhn fc Shapiro! 120071 : lO'Shea Sz Normanl 120081 1 . 
While the Lyman Werner background from early star formation has a primarily negative 
feedback effect, other processes (e.g., cooling i n supernova remna nt shocks) may mitigate 
the suppression of H2 molecules (|Ferraralll998l ; iRicotti et al.l 1200 ll ). 



3. ANALYTICAL CALCULATION 

In this section we describe how we assign the luminosity to the halos and the IGM in 
the simulation. Note that our method is fully analytical, and thus can be adopted to any 
other reionization simulations. 



3.1. Luminosity of the halos 



Luminosity within the halos is dominated by five radiative processes: stellar (black- 
body) emission, and the nebular emission including free-free, free-bound, and two-photon 
emission, as well as any emission lines (here, Lyman-a is the most important one for our study 
of NIRB). The luminosity of each component can be found analytically using the equations 
in FK06. Equations for the stars' luminosity, temperature, number of ionizing photons pe r 
second, and lifetime were based on equations from T able 3 of iFernandez fc Komatsul (120 08). 



which were fit from stellar models or fitting functions ( IMarigo et al.ll200ll ; lLejeune fc Schaerer 
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200ll : ISchaere^ l2002h . 



First, let us define the "volume emissivity." The volume emissivity (luminosity per 
comoving volume per frequency), p{v), is related to the "emission coefficient" (luminosity 



per comoving volume per fr equency per steradian), j u , in ISantos. Bromm &: Kamionkowski 



(120021 ); ICooray et al.l (120041 ) by p(u) = &irj v . In other words, the luminosity is given by 

dL = p(u)dudV = j u dQdudV, (1) 

where dV is the comoving volume element, and dQ is the solid angle element. Integrating 
p(u) over u, one obtains the "comoving luminosity density," pi, as dL = pidV, where 
Pl = J p{y)dv. 

When the main-sequence lifetime of stars under consideration is shorter than the time 
scale at which the star formation takes place, the volume emissivity is given by a product 
of the star formation rate (the stellar mass density formed per unit time), p*(z), and the 
ratio of the mass-weighted average of the total radiative energy (including stellar emission 
and reprocessed light) emitted over the stellar lifetime to the stellar rest-mass energy, (e°) 
(see Eq. (2) of FK06): 

p{v, z) = J>>, z) = p^z)c 2 (2) 

a a 

where 

_ g 2 dm [Lg(m)r(m) / \mc 2 )} f(m)m 
J dmj(m)m 

Here, m is the stellar mass, L°(m) is the time-averaged luminosity per frequency of a given 
radiative process a (which includes the stellar, free-free, free-bound, two-photon, and Lyman- 
a emission), r(m) is the main sequence lifetime, and f(m) is the initial mass spectrum of 
stars under consideration (specified later in § I3.2p . Note that (e°) may also be interpreted 
as a ratio of the total radiative energy within a unit frequency to the total stellar rest-mass 
energy, 

f m2 dmf(m)L?.(m)T(m) 
JJ^ 2 dmf{m)mc 2 

Either way, (e") is a convenient quantity that tells us how much of the stellar rest-mass 
energy is converted into the radiative energy within a unit frequency interval. 

In FK06 we have shown that (e") can be calculated robustly for a given stellar popu- 
lation, i.e., f(m), using the basic stellar physics and radiative processes in the interstellar 
medium. For the radiative processes and stellar populations we consider in this paper, 
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v k e Z) ~ 10 3 ( see Figure 2 of FK06). From this one may obtain a quantity that is commonly 
used in the literature, the luminosity per stellar mass, I", as 

= p a (v, z) = dlnp^z) C dmf(m)L«(m)T(m) 
P*{z) dt f™*dmf(m)m 

In this expression one may identify d In p* (z) / dt as the inverse of the star formation timescale, 
tsp{z), i.e., t$F{z) = [<i In p* (2;) / cit] _1 Ji] Therefore, we finally obtain 

p ^ = 1 C dmf(m)L°(m)r(m) 
tsF(z) f™* dmf(m)m 

when the main sequence lifetime of stars is shorter than the star formation timescale, r(m) < 
tsp(z). In terms of (e") we may also write Eq. (jHJ) as l"(z) = (e")c 2 /^SF(z). 

On the other hand, when the star formation timescale is shorter than the main sequence 
lifetime of stars, £sf(z) < T { m ), we find a different expression for /" (see Eq. (A6) of FK06): 

f m2 dmf(m)L^(m) 

ja J mi J y 1 v\ 1 ,„\ 

J dmj(m)m 

and no longer depends on z as long as f(m) does not depend on z. From Eqs. (jH]) and (J7J) 
we find that the former is roughly r/tsF times the latter. In other words, if one misused the 
latter form when r < tsF, one would over-estimate the signal by a factor of ~ tsp/r, which 
can be as large as a factor of 10 for short-lived, massive stars with ~ 100 M Q E 

For the precise calculation one should use both expressions depending on the situation; 
however, to simplify the analysis, we shall use either Eq. (jBJ) or (J7|), depending on the ratio 
of the stellar lifetime averaged over the initial mass spectrum and weighted by the luminos- 
ity (since more massive, shorter lived stars will contribute more to the overall luminosity), 
(t) = J^ 2 dmf(m)r(m)L/ f™ 2 dmf(m)L, to the star formation timescale (where L is the 



then the 



: assumes that the star formation is triggered by mergers of dark matter halos, 
star formation timescale may be related to the halo merger rate, i.e., tg F (z) = 
[J dMhMh,(d?nh/dMhdt)] /\ f dMhMh(dnh/dMh)], where dnh/dMh i s the mass function of dark matter halos. 
This a pproach was used bv lSantos. Bromm fc Kamionkowskil ( 20021) ; ICoorav et al.l ( 2004 ) ; ICoorav fc Yoshida 
(|2004| ). In this pape r we shall use £sf = 20 Myr as our fiducial value, to be consistent with the value used 



by the simulation of llliev et al.l (|2008bl ). We also study the effects of changing tsp in § 16.11 



Salvaterra fc Ferraral (120031 ) ; iMagliocchetti. Salvaterra. fc Ferraral (120031 ) ; iKashlinskv et al.l (|2004l) used 



Eq. © for r < t SF , 
w t SF /r. 



and thus their predicted amplitudes of NIRB are likely over-estimated by a factor of 
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bolometric luminosity). In the simulations of llliev et al.l (120061 120071 l2008bl ). the star forma- 



tion timescale takes on a universal value, tgp « 20 Myr. For all the stellar populations, the 
luminosity- weighted lifetime is shorter than 20 Myr; thus, Eq. [7] will be our fiducial formula. 

To compute for each radiative process, we use a black-body for the stellar component, 
I* (see Eq. (6) of FK06). This emission is cutoff above 13.6 eV, so all of the ionizing photons 
go into producing emission in the nebula or the IGM. The expressions given in § 2.3, 2.4, 
and 2.5 of FK06 are used for the nebular processes. We then integrate I" over a band of 
observed frequencies V\ to v% to obtain the band- averaged luminosity per stellar mass, l a , as 

r- u 2 (l+z) 

r{z) = / dv r v {z). (8) 
./ 1/1(1+*) 



Following llliev et al.l ( 120061 . 120071 l2008bl ). we assume all halos have a constant mass-to- 
light ratio. With the luminosities per stellar mass, l a {z), computed, we obtain the luminosi- 
ties of the halo, L h (z), by multiplying l(z) by the total stellar mass per halo, f*M h (Q b /Q m ), 
where is the total halo mass (including dark matter and baryons), and /* is the star 
formation efficiency, which is the fraction of baryons that can form into stars over the star 
formation timescale isF- We find 

^ = PC*) + (1 - /esc) [/"(Z) + l f \z) + + ?*"{Z)] } , (9) 

where / esc is the escape fraction of ionizing photons from the halo. Only those photons that 
do not escape into the IGM produce nebular emission within the halo. 

From this result one may conclude immediately that the NIRB power spectrum from 
halos, which is proportional to (L^/M^) 2 , is proportional to f%. Also, Lh/Mh goes down 
as /esc approaches unity, for which all the ionizing photons would escape halos, and thus 
no nebular emission would be left in halos. The stellar properties, such as metallicities and 
initial mass spectra, affect only l a . 



3.2. Stellar Populations 



The simulations from llliev et al.l (l2008bl ) define a quantity, / 7 , which is proportional to 
the number of ionizing photons that escape into the IGM: 



(10) 



where Ni is the number of ionizing photons emitted per stellar atom. When modeling stellar 
populations in our calculations, we shall assure that each of our models agrees with / 7 = 250, 
which was used in the simulations. 
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Population 


Initial Mass Spectrum 


mi,m 2 


(r) (Myr) 


Ni 


fesc 


/* 


Pop III 


Salpeter 


3M , 15OM 


8.08 


5600 


0.22 


0.2 


Pop III 


Larson, m c = 25OM 


3M , 5OOM 


2.45 


25000 


0.1 


0.1 


Pop III 


Salpeter 


3M Q , 15OM 


8.08 


5600 


0.9 


0.05 


Pop III 


Larson, m c = 250M Q 


3M , 5OOM 


2.45 


25000 


1 


0.01 


Pop II 


Salpeter 


3M , 150.1/ . 


9.04 


2600 


0.95 


0.1 


Pop II 


Larson, m c = 5OM 


3M , 150.1/ . 


4.87 


12000 


0.9 


0.023 


Pop II 


Salpeter 


3M , 15OM 


9.04 


2600 


0.19 


0.5 


Pop II 


Larson, m c = 5OM 


3M , 15OM 


4.87 


12000 


0.098 


0.21 



Table 1: Stellar populations ("Pop III" are metal-free, and "Pop II" are metal-poor with the 
metallicity of Z = 1/50 Z & ), parameters for initial mass spectra, the luminosity- weighted 
main sequence lifetime of stars ((t)), the corresponding number of ionizing photons per stellar 
atom (Ni), escape fractions of ionizing photons (fesc), and the star formation efficiency (/*). 
Note that f esc and /* are tuned such that the value of f 1 = fcscf*N is held fixed at / 7 = 250 
which, when combined with the star formation timescale of igp = 20 Myr, can reionize 
the universe such that the resulting electron-scattering optical depth is consistent with the 
WMAP data. 



Iliev et al.l (l2008bl ) have shown that this choice of / 7 , combined with the universal star 
formation timescale of t$F = 20 Myr, can reionize the universe successfully with the resulting 
electron-scattering optical depth consistent with the WMAP data. Within this framework, 
since / 7 /tsF is the only free parameter, models with the same / 7 /isF would produce the 
same reionization histories. (For the simulation case f250C on which our calculations here 
are based, for example, the globally-averaged ionized fraction of the IGM was found to be 
50% at z = 8.3 and 99% at z = 6.6.) To keep / 7 /£sf constant, the star formation efficiency 
must decrease as the escape fraction increases. The various populations that were modeled 
are shown in Table [TJ 

We modeled both zero metallicity stars (Population III; Z = 0) and low metallicity 
stars (Population II; Z = 1/50 Z & ) with either a heavy or a light initial mass spectrum, 
accompanied with either a low escape fraction (/ esc ~ 0.1) or a high escape fraction (/ esc ~ 1). 
While we try to simulate a range of parameters, it is good to keep in mind that our choice of 
/* and / esc for a given / 7 is basically arbitrary. The level of NIRB fluctuations can change 
significantly when paired with different assumptions for the metallicity, mass, and values for 
/esc and /#. 



A lighter mass distribution of stars is represented by a Salpeter initial mass spectrum 
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f lSalpetei]ll955h 



f(m) oc m 



-2.35 



We use mass limits of mi = 3M and m>2 = lbOM^ for th is spectrum. Heavier stars are 
represented by a Larson initial mass spectrum ( lLarsonlll998l ) 



f(m) oc m 1 ( 1 H 



-1.35 



(12) 



with mi = 3M Q , m 2 = 
rri2 = 15OM , and m c 



5OOM , and m c = 25OM for Population III stars and m\ = 3M , 
= 5OM for Population II stars. 



In Figure [T] we show z/Z" (in units of nW Mq 1 ), and in Figure [2] we show l a (z) (in units 
of nW Mq 1 ) averaged over a rectangular bandpass from 1 — 2 fxm, for the stellar populations 
we consider in this paper. In the relevant redshift range, 7 < z < 15, the stellar, two-photon, 
and Lya emission are the most dominant radiation processes, and all of them are on the 
order of I ~ 10 38 nW Mq 1 (20 Myr/t SF ). 



3.3. Luminosity Density from IGM 

Photons that do escape the halos go into producing emission in the HII region surround- 
ing the halo in the IGM (free-free, free-bound, two photon and Lyman-a emission). The 
emission in the HII region can be found using the volume emissivity, p(z/), i.e., luminosity 
per comoving volume per frequency, or luminosity density per frequency (see Eq. ([T|) for the 
precise definition). 

Since all of the radiative processes we discuss in this section are proportional to the 
number density squared, we need to be careful about the comoving versus proper quantities. 
The proper volume emissivity is proportional to the proper number density squared, i.e., 
Pprop oc np rop . As the comoving volume emissivity is p com = a 3 p prop = p pr op/(i + z) 3 and the 
comoving number density is n com = a 3 n prop = n prop / (1 + z) 3 , we obtain p com oc (1 + z) 3 n 2 com . 
This factor of (1 + z) 3 simply reflects the fact that the IGM was denser at higher redshift, 
and thus the IGM was brighter. In the following derivations n always refers to the comoving 
number density. 

For free-free and free-bound emission, the volume emissivity is 

e -hu/kT g 

Pff,fb(^, z) = Atx(\ + z) 3 n e n p ^ c , (13) 

where n e and n p are the comoving number density of electrons and protons respectively, 7 C 
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Fig. 1. — Luminosity spectrum per stellar mass. The stellar, z/Z* (triple-dot dashed red line), 
free-free, vVJ (dotted purple line), free-bound, vl^ (dashed light blue line), two-photon, 
vtp (dot dashed green line), and Lyman-a emission, ul^ ya (solid dark blue line), are shown 
in units of nW Mq 1 as a function of the rest-frame energies. The stars are at z = 10, but 
the redshift affects the profile of t he Ly man-a line only, which was taken from Eq. (15) of 



Santos. Bromm Sz Kamionkowskil (120021 ) . 
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Fig. 2. — Luminosity per stellar mass averaged over a rectangular bandpass from 1 — 2 /im. 
The stellar, I* (triple-dot dashed red line), free-free, tf* (dotted purple line), free-bound, V b 
(dashed light blue line), two-photon, P 7 (dot dashed green line), and Lyman-a emission, 
(solid dark blue line), are shown in units of nW Mq 1 as a function of redshifts. Free- free 
and free-bound both decrease with redshift. This is because both decrease with energy, and 
as redshift is increased, the bandwidth corresponds to higher rest-frame energies. The initial 
rise in Lyman-a is due to the wing of the line. At z ~ 15.5, the line hits the end of the band 
where there is no more Lyman-a emission. Stellar emission increases initially because there 
is more emission from the star as energy increases, and later decreases as the bandwidth 
begins to sample energies above 13.6 eV. Two photon emission is cut off at z ~ 15.5, which 
corresponds to the band sampling above 10.2 eV, above which there is no emission. 
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is the continuum emission coefficient including free-free and free-bound emission: 



Tc 



fk 



9ff 



E 

n=2 



-9fb{n) 



(14) 



where x n = Ry/ (kT g n 2 ), and gbf(n) are the Gaunt factors for free-free (which is thermally 
averaged) and free-bound emission, respectively, fk is the collection of physical constants 
which has a numerical value of 5.44 x 10~ 39 in cgs units, and T g is the gas temperature, 
which we took to be 10 4 K (see § 2.3 of FK06 for more details). 



Using the charge neutrality, n e 



n p , we write 



(15) 



where nn is the number density of hydrogen atoms and X e is the ionization fraction, both 
of which are given in the simulation. The volume emissivity is therefore given by 



4vr(l + z) 3 7c 



-hu/kT g 



1/2 



(16) 



The two-photon emissivity is 

p 2 ~ ( (is, z) = (1 + Z 



(17) 



A fraction of photons that make the 2 — 1 transition, (1 — fhya), go into two photon emission, 
while the remainder, fh ya , produce the Lyman-a line. The precise value of f^a depends 
slight ly on the temperature of gas, and for a gas at 10 4 K the value of f\, ja is 0.64 ( iSpitzer 



19781 1 . Here, «b is the case B hydrogen recombination coefficient given by 

2.06 x 10- 11 



1/2 



-0(T g ) cm 3 s-\ 



where 0(T g ) is given by ISpitzerl (119781 ) . Here, P{y) is the normalized probability per two 
photon decay that one photon is in the range dy = du/u^ ya , which can be fit as (Eq. (22) of 
FK06) 

P{y) = 1.307 - 2.627(y - 0.5) 2 + 2.563(y - 0.5) 4 - 51.69(y - 0.5) 6 , (19) 
and i^Lya is the frequency of Lyman-a photons. Using % and X e , we write the emissivity as 



P2 7 (^, Z) 



; i + ^)3^p (y)(1 _ /LyQ ) aB . 



(20) 
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For Lyman-a, 



Ph ya (v, z) = (1 + z) 3 f Lya hv Lya n e n p a B (j)(is - v Lya ), 



(21) 



where 4>(v — v-\, yn ) is the line pro fi le of the Lyman-a line, given in lLoeb &: Rybickil ( 119991 ) : 
Santos. Bromm fc Kamionkowskil (120021 ) . Using n# and X e , we get 



'1 + Z) fLyahVLyaOtB^iv ~ 



yOL) 



(22) 



p IGM (u,z 



Collecting all the processes we obtain the volume emissivity of the IGM as 

e -hv/kT g 



[I + ZY { 47T 7c 



-,1/2 



a B hv L 



VQ 



(1 — fhya) j L7C ^ + fhya< 



V 



Lya 



(23) 

We are now in a position to find the emission of the IGM by pairing these formulas with the 
hydrogen number densities (tih) and the ionization fractions (X e ) from the simulations. In 
Figure |3] we show vp a (v, z)/[(l + z) 3 n 2 H Xf\ (in units of nW m 3 ) for individual processes as 
a function of the rest-frame energies. 



4. LUMINOSITY-DENSITY POWER SPECTRUM 

The three-dimensional power spectrum of over-luminosity density, <5p^(x), is given by 

(^(k)^ (k')> = (2ir) 3 P L (k)5 3 (k - k'), (24) 

where Pi{k) is the luminosity-density power spectrum, and 5p L (k) is the Fourier transform 
of the over-luminosity density field, 5p^(x). The over-luminosity density field is related to 
the excess in the volume emissivity over the mean, 5p(u,x), integrated over the observed 
bandpass v\ to z/ 2 , as 

/•i/ 2 (l+z) 

5p L (x,z)= / dv Sp(u,x, z). (25) 

Jvx{X+z) 

In the following derivations we do not write z explicitly for clarity. 



4.1. Halo Contribution 

How do we calculate the halo contribution from a given simulation box at a given z7 
For the halo contribution, 5p^ al °, we have 

^r(x) = ( (26) 

\M h J V ce u 
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Fig. 3. — Volume emissivity spectrum of the IGM, up a (i/, z), divided by (1 + z) 3 n 2 H X 2 , for 
individual processes in units of nW m 3 as a function of the rest-frame energies. (Note that 
this quantity does not depend on z.) We use the ionized gas temperature of 10 4 K. Free-free 
(dotted purple line), free-bound (dashed light blue line), two-photon (dot dashed green line), 
and Lyman-a emission (solid dark blue line) are shown. 
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where M ce n(x) is the total mass of halos within a given cell, V ce \\ is the volume of each cell, 
and the bars denote the volume average over the simulation box. Throughout this paper, 
we always include both the stellar contribution as well as the nebular contribution when we 
refer to the "halo contribution." 

Since we assume that halos have a constant mass-to-light ratio, Lh/Mh does not depend 
on x or M h (but it depends on z), and is given by Eq. ([9]). Since we assume a constant mass- 
to-light ratio, the luminosity density <5p^ al ° is linearly proportional to the halo mass density, 
5p\f°, such that (5p^ alo (x) = (L h /M h )5p\f°(x.), where <5p^ alo (x) is the mass over-density of 
halos given by 



r halo / 

°Pm ( 



M celI (x) - M 



cell 



cell 



dM h M h 



dnh{x.) 
dM h 



dn h 
dM h 



(27) 



Here, drih{x)/dMh is the number density of halos per mass within a cell at a location x, and 
dfih/dMh is its average. Therefore, the luminosity-density power spectrum of halos, P^ alo (/c), 
is simply proportional to the mass-density power spectrum of halos, P^f lo (k), as 



jhalo 



(k) 



(28) 



The shape of P£ (k) is determined by that of the halo mass-density power spectrum. In 
other words, one only needs to compute P^f°(k) from simulations, and the analytical calcu- 
lations given in § 13.11 supply Lh/Mh for a given stellar population and observed bandpass. 



Specifically, we compute P^f lo (/c) from the simulation as follows (e.g.. |Jeong fc Komatsu 



2009|): 



(1) Use the Cloud- In-Cell (CIC) mass distribution scheme to calculate the mass density 
field of halos on 256 3 regular grid points, i.e., M ce n(x)/V^ e ii, from the halo catalog. 

(2) Fourier-transform the excess mass density, [ilf cc u(x) — Af ]/V^ e ii, using 

(3) Deconvolve the effect of the CIC pixelization effect. We divide P(k, z) = |c5(k, z)\ 2 at 
each cell by the Fourier transform of the CIC kernel squared: 



w(k) = n 



i=l 



2k N 



(29) 



http: / / www.fftw.org 
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where k = (k\, k 2 , ks), and k^ = tt/H is the Nyquist frequency (H is the physical size 
of the grid). In terms of the number of grids along one axis, iV mes h, one may write 
H = L hox /N mcsh , and 2k N = N mcsh (2iT / L hox ) = N mcsh Ak, where A A; = 2%/L hox is the 
fundamental frequency of the box. (For our simulation, L^ ox = 100 /i -1 Mpc.) We also 
try a different deconvolution scheme that attempts to reduce the aliasing effect (IJing 



2005|): 



3 r 



W(k) = ] 



sm 



2k 



N 



(30) 



We then use Pm(^) up to /c max 
the same answer. We find k m ~ x 



below which both of the deconvolution schemes yield 
- 5 Mpc" 1 . 



(4) Compute Pj^ik^z) by taking the angular average of CIC-corrected P(k,z) within a 
spherical shell defined by k - Ak/2 < |k| < k + Ak/2. 



In the previous work on NIRB fluctuations (IKashlinsky et al.ll2004l ; ICooray et al.ll2004l ) 
the linear bias model was used, i.e., P£f lo (/c) was assumed to be linearly proportional to 
the underlying (linear) matter power spectrum. However, for such high redshifts halos are 
expected to be highly biased, and thus non-linear bias cannot be ignored. In other words, it 
is no longer correct to assume that P\f°(k) is linearly proportional to the underlying matter 
power spectrum. 

To study this further, in Figure|4]we show P^f°{k) (in units of Mj, Mpc -3 ). Also shown 
in Figure|4]is the shot noise, P]$ ot , where P^j ot = / dMhM^dfih/ dMh (dfih/dMh is the mean 
halo mass function), the linear matter density fluctuations times the mean mass density 
squared (Pn n (k)(p^f°) 2 ), where p^ 10 is the mean mass density of halos within the simulation 
box, and the bias, given by: 



bes(k) 



P M°( k ) ~ Pii(k) 

i^ypUk) 



(31) 



By comparing P\f°{k) (with the shot noise subtracted) with the power spectrum of 
linear matter density fluctuations times (p^' ) 2 ' we fi 11 ^ that, on large scales (A; < 0.1 Mpc -1 ), 



they are related by P M 



ia,lo (k)/(p\f°) 2 fa b\P^ n (k) with the linear bias factor being b\ ~ 5 at 
10, a highly biased population. The bias increases monotonically as 



z = 6 to bi ~ 10 at z 

we go to smaller scales, significantly boosting the power in the halo distribution relative to the 
matter distribution. This changes the prediction for the shape of the ang ular power spectrum 
qualitatively, compared with the previous results given in the literature ( ICooray et al.ll2004l : 



Kashlinsky et al.ll2004l ). This behavior of non- linear bias with redshift is consistent with that 
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Fig. 4. — Non-linear bias of the halo mass-density power spectrum. (This is not the 
luminosity-density power spectrum; see § 14.11 for the precise definition.) Top left panel: 
The power spectra of the halo mass density, P\f°(k) are shown as the solid lines (z — 6 to 
10 from top to bottom), the linear matter power spectra times the mean halo mass density 
squared, Pn n (k)(p\f°) 2 , are the dashed lines, and the shot noise power spectra, P^ ot , are 
the dotted lines. Top right panel: We show the bias, ^/[P^f°(k) - P$ ot (k)]/[(p\f°) 2 P rm (k)], 
(z = 10 to 6 from top to bottom). The bias increases significantly as we go to smaller scales, 
and this effect has been ignored in the previous calculations of the power spectrum of NIRB 
fluctuations. Note that the minimum halo mass resolved in the simulation is 2.2 x 10 9 M & . 
The de gree of non- linear b ias would be smaller for a smaller minimum mass (see, e.g., Fig- 
ure 6 of lTrac fc Cenl 120071 ). Bottom left panel: The linear power spectrum, k 3 Pu n {k) / (2ir 2 ) . 
Bottom right panel: Same as top right panel, but on a log-log axis. 
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5 10 15 20 



Fig. 5. — 5 c /a(MMin, z) versus redshift. The halos resolved in our simulation, with M > 
M m[Q = 2.2 x 10 9 Mq, are located on rare peaks (S c / a(M min , z) > 2.5) at z > 7. 



expected from the halo model (jCooray fc Shethll2002l ). These halos are very rare, located on 



high peaks with 5 c /a(M min , z) > 2.5 (see Figure [5]). 
This motivates our writing P^ alo (k) as 



2 



phalo^ = blMP^k), (32) 

where the pre-factor, p\f°Lh/Mh, is the mean halo luminosity density. In the left panel of 
Figure [6] we show p h ^°L h /M h (in units of nW Mpc~ 3 ) function of redshifts. We find 
that the redshift evolution of p^J^Lh/Mh is very rapid; thus, the redshift evolution of the 
halo luminosity density power spectrum, P^ alo (k), is dominated by that of the mean halo 
luminosity density. 



What determines the evolution of the mean halo luminosity density? The answer is 
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Fig. 6.— {Left) Mean halo luminosity density computed from our simulation, 
p\f°(z)Lf l (z)/Mf l , where l"(z) is from equation in units of nW Mpc -3 as a function 
of redshifts, for various stellar populations given in Table [U The waves in the lines where 
/esc are higher are from the discrete redshift sampling of the Lyman-a line. We averaged 
the luminosity over a rectangular bandpass of 1 — 2 /im. (Right) Halo mass collapse frac- 
tion, pjf°(z)/(Q m p c o), as a function of redshifts. The redshift evolution of p\f°Lh/Mh is 
essentially determined by that of p^f 10 . 

simple: it is determined by the rate at which the mass in the universe collapses into halos. 
To show this, in the right panel of Figure [6] we show the halo mass collapse fraction, or 
the ratio of p\f° to the mean comoving mass density of the universe, Q m p c0 , where p c0 = 
2.775 x 10 11 h 2 M Mpc -3 is the critical density of the universe at the present epoch. The 
evolution of the collapse fraction is fast, explaining the fast evolution of the mean halo 
luminosity density. 

As halos are discrete objects, and we do not expect to resolve individual halos con- 
tributing to the diffuse NIRB, the observed NIRB power spectrum is a sum of the clustering 
component and the shot noise component. If the shot noise dominates over the clustering 
component, it would be very difficult to ascertain information on the structure from the sig- 
nal of the NIRB. The shot noise component can be estimated by integrating the luminosity 
squared over the mass function: 

pr = (£)V = / iM h M>*, (33) 

where we have again assumed that each halo has a constant mass-to-light ratio, i.e., L h /M h 
is independent of Mh- 
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4.2. IGM Contribution 



For the IGM contribution, we have 



5p 



IGM, 



P 



IGM 



n 2 ^ 



C ce ii(x)r^ cll (x)X e 2 u (x) - (C ce un 2 ell X 2 u ) 



(34) 



where p 
i.e., p IGM 



IGM 



is the volume emissivity of the IGM, integrated over the observed frequencies, 



dvjr (y), C ce n, ^ceib ^e,ceii are the clumping factor, the comoving 
number density of hydrogen atoms, and the ionization fraction within a cell, respectively. 
We compute n ce ii using 

^6 PM.cell 



- "o rju, ceil / or \ 
"-cell = 7^ , (35) 

\l m firrip 

where \x = 0.59 and m v are the mean molecular weight of ionized gas and the proton mass, 
respectively. We have used the mass density of iV-body particles per cell, pM.ceii, multiplied 
by the baryon fraction, Qf,/Q m , for computing the mass density of baryons per cell, as we 
have assumed that gas traces dark matter particles, i.e., iV-body particles. The clumping 
factor, C ce ii = n 2 actual /n 2 eU , relates the actual density squared to the square of the density 
averaged within a cell. In other words, C ce ii captures the sub-grid clumping that is not 
resolved by the simulation. 



Following llliev et al.l (120071 ). we make a simplifying assumption that C cc ii takes on 
the same value everywhere in the simulation, and evolves with redshift z as C ce \\(z) = 
26.2917e- ai8222+0 - 00350522 ; thus, we have 



8p 



IGM, 



26.2917e 



-0. 1822^+0. 0035052 2 / P 



JGM 



n 



cell 



x)X 



e.ccll \ 



in 



cell e,celU 



(36) 



Note that p M /{n 2 H Xl) does not depend on x, and is given by Eq. fl23|) integrated over a 
rectangular bandpass of 1 — 2 fim in the observer's frame. 



5. RESULTS 

5.1. Luminosity-density Power Spectrum 

In Figures [7] and |S] we show the luminosity-density power spectra, Pi{k), for halos and 
their associated HII regions in the IGM for two of our populations: Population II stars with 
a Salpeter initial mass spectrum with / osc = 0.19 and /* = 0.5 (Figure |7J) and Population III 
stars with a Larson initial mass spectrum with / esc = 1 and /* = 0.01 (Figured]), assuming 
a rectangular bandpass from 1 — 2 /xm. 
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Fig. 7. — (Left) Luminosity-density power spectrum of halos with Pop II stars obeying a 
Salpeter initial mass spectrum, / esc = 0.19, and /* = 0.5, assuming a rectangular bandpass 
from 1 — 2 /im. The shot noise for the halo contribution is also shown as the dotted lines. 
(Right) Luminosity-density power spectrum of the IGM. The ionization fraction of the IGM 
reaches 0.5 at about z ~ 8.3. On large scales where the shot noise is sub-dominant, we find 
P L (k) oc k~ 3/2 , which yields Q oc Z~ 3 / 2 or l 2 Q oc I 1 ' 2 (see §E2). 
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Fig. 8. — (Left) The same as the left panel of Figure H with Pop III stars with the Larson 
initial mass spectrum, / csc = 1, and /* = 0.01. (Right) The same as the right panel of Figure 
[7] for comparison. 
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The luminosity-density power spectra of halos are approximately power-laws over the 
entire range of wavenumbers that the simulation covers. At the highest redshift bin, z ~ 16, 
the power spectrum is entirely dominated by the shot noise at all scales. The lower the 
redshifts are, the more power in excess of the shot noise we observe on large scales (because 
the shot noise is most important on small scales). The growth of the power spectrum is 
partly driven by the growth of linear matter fluctuations as well as that of halo bias, i.e., 
the clustering of halos is biased relative to the underlying matter distribution. As we have 
shown in the previous section, the bias of halos that we observe in the simulation is highly 
non-linear, and thus has an important implication for the predicted shape of the observed 
power spectrum of NIRB fluctuations. However, as we have shown in § 14.11 the evolution 
of Pi{k) is almost entirely driven by the fast growth of the mean halo luminosity density, 
p\f°(z)L h (z)/M h (see Eq. (132]) and the left panel of Figure [6]). As a result -Pl(^) grows by 
about six orders of magnitude at k — 0.1 Mpc" 1 from z ~ 13 to z ~ 6, which is much faster 
than the growth expected from the growth of bias times the matter power spectrum. 

The luminosity-density power spectrum of the IGM increases quickly as the mean ioniza- 
tion fraction, X e , approaches 0.5 (at about z ~ 8.3 for this particular simulation), especially 
on larger scales. As the ionization fraction increases, the luminosity of the HII region would 
also increase (because luminosity is proportional to X 2 ). Moreover, since we are looking 
at the over-luminosity- density power spectrum of the IGM, the greatest power results when 
there is the greatest difference between luminous regions and the average luminosity of the 
IGM; thus, the power spectrum of X 2 n 2 grows rapidly as X e approaches 1/2. However, this 
rapid growth of the power stops when the entire IGM is ionized (X e = 1), in which case the 
over- luminosity-density power spectrum of the IGM is simply proportional to n 2 . 

The most interesting feature of the luminosity- density power spectrum of the IGM is a 
"knee" feature, which is at k ~ 2 Mpc -1 at z ~ 16, and moves to k < 1 Mpc" 1 at z < 10. 
This "knee" is caused by the typical size of HII bubbles: the knee wavenumber is inversely 
proportional to the typical size of the bubbles. At the highest redshift bin, z ~ 16, the 
bubbles are nearly Poisson-distributed, and thus the power spectrum is flat up to the knee 
scale, k ~ 2 Mpc" 1 , beyond which the power decreases as one is looking at the scales inside 
the bubbles, which are smooth. As the redshift decreases, the knee scale moves to larger 
scales, signifying a growth in the ionized bubbles with time until they merge. At the same 
time, the large-scale power also grows, and the shape of the HII region power spectrum is 
basically the same as that of the halo power spectrum, as the bubbles are created around 
the halos. 



Note that llliev et all ( 12006k 120071 ) studied the power spectrum of ionized gas density, 
and observed a similar trend. The power spectrum of the luminosity density that we have 



-24- 



presented here is the four-point function of the ionized gas density (as the volume emissivity 
is proportional to the ionized gas density squared), and thus it is different from the power 
spectrum of the ionized gas density (which is quadratic in density). 



5.2. Angular Power Spectrum of NIRB Fluctuations 

What about the observable, the angular power spectrum of NIRB fluctuations, Cp. 
We compute the angular power spectrum of NIRB fluctuations, Ci, by projecting Pi(k) on 
the sky. We do this using Limber's approximation, and obtain (see Appendix |A] for the 
derivation) 

c f dz ( I \ 

Cl = j^fJ H (z)r*(z)(i + z y Ph y = V(z)' z ) 1 (37) 

where r(z) = c f Q dz'/H(z') is the comoving distance. We integrate Eq. (13?]) over the range 
of redshifts that the simulation covers for both halos and the IGM, z = 6.0 — 15.7. 

In Figure [9] we show 1(1 + l)Ci/(2ir) for halos with Population II stars with a Salpeter 
mass spectrum and /* = 0.5 (the angular power spectrum for halos with the highest ampli- 
tude) and for Population III stars with a Larson mass spectrum and /* = 0.01 (the angular 
power spectrum for halos with the lowest amplitude), along with the angular power spec- 
trum of the IGM. The halo contribution at small scales, i.e., I > 10 4 , is comparable to the 
shot noise contribution. When the shot noise is subtracted (see Figure [TO}) , we find that 
1(1 + l)Ci/(2ir) is nearly a power-law, 1(1 + l)C;/(27r) oc /°' 5 , with no sign of a turn-over, 
which would be expected from the shape of the project ed linear matter power spectrum. This 



is in a stark contrast with the previous calculations (IKashlinsky et al.l 120041 ; ICooray et al. 



20041 ). which predicted a turn-over at / ~ 10 3 . They assumed that the luminosity-density 
power spectrum was given by the linear bias model, in which the halo power spectrum is a 
constant times the matter power spectrum. Our calculations, which are based on a realistic 
simulation, indicate that the simple linear bias model is not valid for these populations. This 
is expected, as these populations are very highly biased, and therefore non-linear bias must 
also be large, as demonstrated already in Figure HI 

On the other hand, there is no freedom in changing the amplitude of the IGM power 
spectrum for a given simulation, i.e., a given / 7 /isF5 thus, we show only one line for the 
IGM contribution in Figure M (the lowest line). For the parameter space explored here, 
the halo contribution can be as low as being only slightly over an order of magnitude (for 
PopIII Larson with / esc = 1 and /* = 0.01) to about 10 6 times greater (for PopII Salpeter 
with / esc = 0.19 and /* = 0.5) than the IGM contribution. If we were to increase / 7 
(which is possible using additional simulations in future work, although one has to make 
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Fig. 9. — Angular power spectra of NIRB fluctuations, Ci, from halos in comparison to the 
angular power spectrum of the IGM (the bottom line). We show C\ from halos that have 
Population II stars with a Salpeter mass spectrum and /* = 0.5 (the angular power spectrum 
that has the highest amplitude) and Population III stars with a Larson mass function and 
/* = 0.01 (the angular power spectrum with the lowest amplitude and which is closest to the 
angular power spectrum of the IGM). The dotted lines show the level of the shot noise. In 
Figure [LTJ we show how the amplitude of the power spectrum changes between populations 
with various escape fractions of the ionizing photons into the IGM, / esc , and star formation 
efficiencies, /*. (Right panel) Same as the left panel, except divided by The IGM 
contribution is not shown. 
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Fig. 10. — The angular power spectrum from the clustering of halos (solid line), i.e., the 
angular power spectrum minus the shot noise contribution. The dotted line has a slope of 
/°' 5 . The clustered angular power spectrum shows no evidence of a turnover that was claimed 
to exist in the literature. This is because previous analytical models in the literature based 
their power spectrum on the linear bias model, which is not valid for this population, which 
has a high level of non-linear bias. The minimum halo mass used in this calculation is 
2.2 x 10 9 Mps. 
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Fig. 11. — (Left panel) The change of the angular power spectrum as a function of the escape 
fraction, / esc , for our selected samples of stellar populations. Each amplitude is scaled with 
relation to the angular power spectrum of Population II stars with a Salpeter mass spectrum 
and /* = 0.5. (Since the shape of the angular power spectra are the same for all stellar 
populations, this ratio is the same for all wave numbers.) (Right panel) The dependence of 
the angular power spectrum on /*. The solid line shows C\ oc fl. Note that each stellar 
population has a different set of /*, f csc , and N iy and thus both panels show a slice of the 
multi-parameter space. 

sure that the resulting electron-scattering optical depth is consistent with the WMAP data), 
a wider range of parameters / esc , /* and iV, could result. This is a good news, as this gives 
us an opportunity to study the physics of the reionization using the power spectrum of 
NIRB fluctuations. Sensitive surveys may be able to detect a change in the shape of the 
power spectra that would be a result of the IGM power spectrum. This may be one way of 
constraining / esc observationally. 

In Figure [TTJ, we show the amplitude of the angular power spectra of other stellar 
populations scaled to the angular power spectrum of Population II stars with a Salpeter 
mass spectrum and /* = 0.5. As we find in Eq. Qj, the luminosity-density power spectrum 
of halos is about proportional to f%, and one of the terms in the power spectrum (nebular 
contribution; the second term in Eq. Q) depends on (1 — / csc ). Therefore, for a fixed 
f-y — fcscf*Ni and fixed iVj (i.e., fixed stellar population), the angular power spectrum of 
the halo contribution must always increase as we increase /*, as increasing must be 
accompanied by the corresponding reduction in / esc , both of which will increase the power 
spectrum of the halo contribution. 

As Ci oc f%, the parameter combinations that maximize /* tend to give the largest 
C[. For a fixed / 7 = f eB cf*Ni this means a lower N iy i.e., lighter mass spectra with larger 
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metallicity (see the 5th column of Tabled]), and a lower / esc . In reality, however, we should 
also take into account the fact that heavier mass spectra produce more luminosity per stellar 
mass, i.e., more Tin Eq. OH]). These factors explain the dependence of the predicted amplitudes 
of l(l+l)Ci/(2n) (averaged over A = 1—2 /xm) on parameters shown in Figure [TTl Populations 
with higher / e sc have lower angular power spectrum. This is to be expected, because as / esc 
increases, less photons are available to create luminosity within the halo. 



6. VARYING THE MODEL: HALO CONTRIBUTION 

In this section, we will focus on the halo contribution to the angular power spectrum of 
NIRB fluctuations, and explore the effects of changing various parameters. 



6.1. THE EFFECT OF THE STAR FORMATION TIMESCALE 



As mentioned in section I3.1[ the star formation timescale will affect the amplitude of 
the angular power spectrum. We have assumed in this work a constant star formation 
timescale of tsF = 20 Myr to make a consistent comparison between the halo and the IGM 
contributions. However, the amplitude of C% from halos depends sensitively on this rather 
uncertain timescale, as the luminosity of halos is proportional to igp, and thus Ci oc l/t| F - 
Motivated by this, in this section we consider two other possibilities: 1) The star formation 
time scale is shorter than the lifetime of the stars, in which case we will use equation [7] to 
compute the luminosity per mass, and 2) the star formation is triggered by mergers, i.e., 



/ dM h M h (d 2 n h /dM h dt) 
J dM h M h (dn h /dM h ) ' 



(38) 



where drih/dMh is the mass function of dark matter halos. For the Press-Schechter mass 
function, we can calculate £sf(^) analytically from 



^SF 



H(z) 



dlnD 



cHn(l + z) 



81 



D 2 (z)a 2 (M n 



H(z)n 



0.55/ 



D 2 (z)a 2 (M n 



(39) 

where 5 C = 1.68, D(z) is the growth factor of linear matter density fluctuations normalized 
such that D(0) = 1, a(M min ) is the present-day r.m.s. matter density fluctuation smoothed 
over a top-hat filter that corresponds to the minimum mass M min , and Q m (z) is the matter 
density parameter at a given z. Note that interpreting this quantity as a merger timescale 
makes sense only when we study the density peaks above the r.m.s., i.e., 5 c /[D(z)a(M m i Q )] > 
1. (Otherwise tsF becomes negative.) This formula has a clear physical interpretation: for 
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a density peak of order the r.m.s. mass density fluctuation, 5 c /[D(z)a(M m i n )] -lwl, the 
merger timescale is of order the Hubble time, i.e., tsF ~ H~ l (z). The higher the peaks are, 
the shorter the merger timescale becomes; thus, in this model, high-z objects (for a given 
mass) have shorter star formation timescales, and are brighter. 

As the reionization history depends on / 7 /isF, changing only tsF without the corre- 
sponding change in / 7 results in a different reionization history. For example, increasing tgp 
by a factor of 10 makes individual sources fainter by a factor of 10, and thus it would result 
in a much slower reionization history. To compensate this one would have to increase / 7 
by a factor of 10. Moreover, if we reduce tsF by a large factor, it would make individual 
sources brighter by a large factor, to the point where we might start detecting these sources 



individually, e.g., as Lyman-a emitters ( [Fernandez fc Komatsul 120081 ) 



In this section, however, we explore the effects of tgp for a given / 7 , to show how 
important this quantity is for predicting the amplitude of NIRB fluctuations without any 
extra information on reionization from WMAP or Lyman-a emitters. 

The angular power spectrum for various assumptions for the star formation timescale 
is given in Figure [T2j The angular power spectrum with the highest amplitude corresponds 
to when the star formation time scale is shorter than the lifetime of the stars. If the star 
formation timescale is given by the merger time of halos (Eq. [39|) . the star formation timescale 
varies with redshift and we obtain the lowest amplitude for the angular power spectrum, as 
the merger timescale at a given redshift is usually comparable to the age of the Universe at 
the same redshift. Our assumption of tsF = 20 Myr lies between these two extremes. 

We can further quantify the uncertainty in C\ from tgp by looking at the mass-to-light 
ratio of the galaxies (see Figure [T^]) . We know very little about the nature of high-z galaxies 
contributing to NIRB. We don't know what the mass-to-light ratio is for these populations. 
An uncertainty of a factor of 100 in the star formation timescale will correspond directly to 
an uncertainty in the mass-to-light ratio of 100, and an uncertainty of 10 4 in the angular 
power spectrum. Early galaxies could be starbursts, with a mass-to-light ratio of less than 
0.1 to 1, or normal galaxies, with a mass-to-light ratio of > 10. The amplitude of C/ is, 
among other things, a sensitive probe of the nature of high- z galaxies. 



6.2. THE EFFECT OF CHANGING z cnd 

The angular power spectra will also depend on what we choose for the end of the star 
formation epoch, z end . The effect of our choice of Zbcgin is minimal, because at high redshift, 
the halos are smaller and dimmer, contributing less to the angular power spectrum. (See 
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Fig. 12. — The effect of the star formation time scale on the angular power spectrum. We 
find the largest amplitudes when tsF is shorter than the main sequence lifetime of stars, 
whereas we find the lowest amplitudes when tsF is given by the timescale of halo mergers 
(Eq. 1551) . The uncertainty due to the star formation time scale is large and can lead to an 
uncertainty in the angular power spectrum of a factor of ~ 10 4 . This reflects our uncertainty 
in the mass to light ratio of galaxies that contribute to the NIRB. However, note that not 
all scenarios shown here yield the reionization histories that are consistent with the WMAP 
data and the abundance of Lyman-a emitters. (Right panel) Same as the left panel, except 
divided by /jf. 



Bolometric Mass to Light Ratio 



Bolometrlc Mass to Light Ratio x f. 



2 1.000 r 



Salpeter, Pop2, f.=0.5 
Larson, Pop 5, f.=0.01 



Fig. 13. — The bolometric mass-to-light ratio for halos for various star formation timescales. 
Uncertainty in the amplitude of the star formation time scale can be equated to the un- 
certainty in the mass-to-light ratio, i.e., the nature of high- z galaxies contributing to the 
NIRB. The upper and lower sets of lines show the PopIII Larson and the PopII Salpeter, 
respectively. (Right panel) Same as the left panel, except multiplied by /*. 
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Figures [7] and [HI) Since halos and IGM will contribute more to fluctuations at lower redshifts, 
we find that the angular power spectrum dramatically drop as we stop star formation at 
higher redshifts (see Figured!]). 

The shape of the angular power spectrum also changes as we vary z en d- As z en d increases, 
the angular power spectrum of the halos steepens. The shape of the angular power spectrum 
from the IGM can also affect the overall slope of the observed angular power spectrum if the 
halo contribution is close to that of the IGM contribution. If the escape fraction is small, 
this effect in the change of shape from the IGM will be less than if the escape fraction is 
large. When z en d is very large, the amplitude of the angular power spectrum of the IGM 
could even be higher than that of the halos. 



6.3. LYMAN-« ATTENUATION 

The Lyman-a line can be attenuated by dust or neutral hydrogen. To understand this 
effect one would have to perform detailed calculations of the radiation transport of Lyman-a 
photons, including scattering of Lyman-a photons; however, such calculations are usually 
quite complex and time-consuming. Therefore, in this subsection we study the extreme 
limit of attenuation: the case where all of the Lyman-a photons are absorbed or extinct. 
How would this affect the angular power spectrum? The effect of the complete Lyman-a 
attenuation is shown in Table [2j The effect of the Lyman-a attenuation is the greatest when 
the Lyman-a line is the strongest (for heavy Pop III stars) and when the escape fraction 
is smaller (so more photons stay within the halo to produce nebular emission). The effect 
of Lyman-a attenuation in the IGM is the highest, because normally a higher fraction of 
emission is coming from the Lyman-a line (in the halos, there is also stellar emission). 



7. COMPARISON TO PREVIOUS WORK 



Cooray et al.1 (120041 ) made fully analytic predictions of the angular power spectrum 
in the NIRB luminosity expected from the first stars in halos. They ignored the IGM 
contribution, which we found to be small relative to the halo contribution for a range of 
parameters we have explored in this paper. They modeled halos with 300 solar mass stars 
for two cases: (1) an optimistic scenario - star formation in halos above 10 5 K, halos forming 
stars from z = 10 — 30, and a star formation efficiency of 100%; and (2) a pessimistic 
scenario - star formation beginning at 5000 K (so the bias is lower), halos forming stars from 
z = 15 — 30, and a star formation efficiency of 10%. Using the same stellar masses (300 M ), 
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Fig. 14. — The angular power spectrum for halos and IGM as z e nd is varied. We show 
the angular power spectrum for the halos with the highest and lowest amplitude of the 
angular power spectrum (Population II stars with a Salpeter mass spectrum and /* = 0.5 
and Population III stars with a Larson mass spectrum and /* = 0.01 respectively) and the 
IGM. The angular power spectrum as shown throughout the rest of the paper has z en d = 6. 
As z cn d increases, the angular power spectrum drops. At very high redshifts, the angular 
power spectrum of the IGM is higher than some of the angular power spectrum of the halos. 
(Right panel) Same as the left panel, except divided by The IGM contribution is not 
shown. 
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Population 


Initial Mass Spectrum 


fesc 


/* 


Cj, Lya atten/ Cz, no atten 


Pop III 


Salpeter 


0.22 


0.2 


0.848 


Pop III 


Larson 


0.1 


0.1 


0.632 


Pop III 


Salpeter 


0.9 


0.05 


0.975 


Pop III 


Larson 


1 


0.01 


1 


Pop II 


Salpeter 


0.95 


0.1 


0.995 


Pop II 


Larson 


0.9 


0.023 


0.974 


Pop II 


Salpeter 


0.19 


0.5 


0.926 


Pop II 


Larson 


0.098 


0.21 


0.825 


IGM 








0.448 



Table 2: The effect of Lyman-a attenuation on the angular power spectrum. Here, we assume 
complete attenuation (no production of Lyman-a photons). The angular power spectrum is 
only slightly affected in most cases, and is more affected in cases where the Lyman-a line 
was strong to begin with (such as heavy Pop III stars). The effect of Lyman-a attenuation 
in the IGM is the highest, as the IGM does not have the stellar contribution, and is mainly 
dominated by the Lyman-a and two-photon emission. 



we ha ve compared our results from the simulation to the optimistic case from ICooray et al. 
( 120041 ) for two different esc ape fractions, and 1, and show the results in Figure [T5] for 
different wavelengths. As in lCooray et al.l ( 120041 ). we use the star formation time scale given 
by the merger time scale (see Eq. 138]) . The angular power spectrum here is 



cr 



dz 



(4vr) 2 J H(z)r 2 (z)(l + z) 2 



Pp ^(1 + Z), V\\ + Z)-k= -ly, Z^j 



(40) 



which gives the angular power spectrum at only one wavelength (rather than that averaged 
over a certain bandpass). The difference between this equation and Eq. ( )37|) is a factor of 
(1 + z) 2 since we are no longer integrating over a range of frequencies (see Appendix lAl for 
the derivation). Note that we do not show C\ at 1 /im: at 1 /im, the emission comes from 
photons that are more energetic than hv = 13.6 eV in the rest frame at z > 10. Because of 
this, there should be no emission from the halos themselves, if one considers halos at z > 10. 
(There would be contributions if one considered halos at lower redshifts, say, z > 6.) 

Since there were not enough halos in our simulation to create an accurate power spec- 
trum a bove z = 16 - 6, our population of stars only goes from 10 < z < 16.6, while the model 



from ICooray et al.l ( 120041 ) included star formation from 10 < z < 30. However, this should 
not make too much of a difference, because halos at higher redshift do not contribute as much 
to the angular power spectrum. In Figure [15] we show the angular power spectrum minus 
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Fig. 15. — Comparison to ICooray et al.l ( 12004] ) (shown as triple-dot dashed lines). We show 
1(1 + l)C//(27r) where C x = v 2 C\ v (see Eq. [40]), from halos at z > 10 that host only very 
massive stars with 300 M , at various wavelengths. The total angular power spectra from 
this work are shown as solid lines, shot noise is shown as dotted lines, and the clustered 
angular power spectra, which are the total power minus the shot noise components, are 
shown as dashed lines. Note that the amplitudes of C\ shown here are much smaller than 
those shown in the previous figures (despite a high star formation efficiency, /* = 1), as we 
have removed t he most dominant, lo wer redshift contributions, z < 10, in this figure, to be 
compatible with lCooray et al.l (120041 ) . See Figure [T41 for the effects of changing the minimum 
redshift of star formation. The mean intensity, for this population of stars at 2/im and 
4/im are 63 and 16 nW m~ 2 sr" 1 respectively, which is already ruled out by observations 
(see section [9]). 
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the shot noise, which will give us the angular pow er spectrum o f the c lustered component, 
which is directly comparable to the quantity from ICooray et al.l (120041 ) . We have included 
all the nebular processes including the free-bound and two-photon emissio n, which are im- 
portant to the overall luminosity of the halo and which ICooray et al.l (120041 ) h ave neglected 



The o verall amplitude of our angular power spectrum is lower than that which ICooray et al. 
( 2004 ) predicted, by a large factor, 10 3 . 

In addition, the angular power spectru m from ICooray et al.l (120041 ) peaks at about 
I ~ 1000 and then turns over. This is because ICooray et al.l (120041 ) did not take into account 
the nonlinear bias in the halo power spectrum. Nonlinear bias will increase the power at 
small scales, especially at high redshifts, where galaxies were more highly biased. We again 
refer to Figure HI which shows the importance of non-linear bias. This greatly affects both 
the amplitude and the shape of the angular power spectrum of the NIRB and should be 
included. 



8. OBSERVING THE FLUCTUATIONS IN THE NEAR INFRARED 

BACKGROUND 



Interpretation of the NIRB data can be a challenging task. Instrument emission, fore- 
grounds and zodiacal light must all be taken into account. Foreground stars and low-redshift 
galaxies, in addition to very faint and the dim wings of galaxies, must be removed. Much of 
the differences in the existing measurements of the fluctuations from stars at high redshift 
result from differences in how lower redshift galaxies are accounted for. Foreground galax- 
ies are removed down to a limiting magnitude (which is usually different between different 
observations). Galaxies fainter than this are taken into account using different methods. 



There have been several observations of the NIRB. lKashlinsky fc Odenwaldl ( 120001 ) found 
fluctuations at the wavelengths from 1.25 to 4.9yum in the images taken by the Diffuse Infrared 
Background Experiment (DIRBE) on Cosmic Background Explorer (COBE) , which were not 
consistent with the Galactic emission or instrument noise. iMatsumoto et al.l ( 120051 ) observed 
the NIRB using the Infrared Telescope in Space (IRTS). They detected a clustering excess 
on scales of about 100' from 1.4 to 4 /im, and an indication of a spectral jump from the 
high redshift Lyman cutoff. This jump could indicate that Population III star formation 



4 This difference may be explained by the fact that lCoorav et al] (|2004l ) actually re scaled the overall ampli- 
tude to fit the mean intensity measured by the Infrared Telesc ope in Space (IRTS) (Matsumoto et al.ll2005l) 



and the Diffuse Infrared Background Experiment (DIRBE) (|Kashlinskv fc Odenwaldl I2000T ) . (A. Cooray, 
private communication.) 
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ended at about a redshift of z ~ 9. Excess fluctuati ons were detected , possib l y from high 
redshift galaxies, at about 1/4 of the mean intensity. iKashlinsky et al.l (|2007d . 120051 ) made 
observations of the fluctuations of the NIRB using the Infrared Array Camera (IRAC) on 
the Spitzer Space Telescope at 3.6, 4.5, 5.8 and 8 /im. Sources were removed by clipping 
pixels containing > 4<r peaks, as well as removing fainter sources identified by SExtractor 
and convolved with the appropriate point spread function of IRAC. Since zodiacal light is 
not fixed in celestial coordinates, it was removed by taking observations six months apart in 
fields rotated by 180°. They detected excess fluctuations (0.1 nW m -2 sr _1 at 3.6 /im) that 
were not consistent with instrument noise, dim wings of galaxies, zodiacal light, or galactic 
cirrus. They claim that it is possible that the excess fluctuations came from high redshift 
galaxies (z > 6.5) or faint, low redshift galaxies. However, since these fluctuations show 
little (< 10~ 3 ) correlation with the ACS source catalog maps, and the power spectrum of 
fluctuations is inconsistent with the Hubble Space Telescope Advanced Camera for Surveys 



(ACS) ca talog galaxies, they stat e it is unlike 
galaxies ( IKashlinsky et al.ll2007al). However, 



of the fluctuations detected by IKashlinsky et al. 



y that these fluc t uations are from faint, low- z 



Thompson et al.l (l2007bl ) claim that the color 



( 12007d . 120051 ) are consistent with objects at 



z < 10, and therefore not from a population of high redshift stars. 



Cooray et al.l (120071 ) observed the NIRB using IRAC at 3.6 /im. They masked the image 
to cut out faint, low redshift galaxies. In their most extensive masked image, they masked 
IRAC sources down to a magnitude of 20.2 in addition to galaxies in ACS catalog. They 
also discarded pixels that had a flux 4<r above the mean. 



Thompson et al.l (j2007al ) also made observations of the fluctuations of the NIRB using 
the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) camera on the Hubble 
Space Telescope at 1.1 and 1.6 /im. The effects of zodiacal light were removed by dithering 
the camera. After removing galaxies down to the fainter ACS and NICMOS detection limit, 
fluctuation power dropp ed two orde rs of magnitude in comp arison to an earlier paper by 
Kashlinsky et al.l (120021 ) . Therefore. iThompson et al.l (j2007al ) confirmed that the observed 
fluctuations reported by IKashlinsky et al.l ( 120021 ) in the 2MASS data are from low redshift 
galaxies (z < 8) (although they are unable to rule out contributions from galaxies in 8 < 
z < 13). Yet, they concluded that an excess fluctuation power in the NIRB of about 
1 — 2 nW m -2 sr _1 could still be from the first stars. Their methodology would miss 
fluctuations that are flat on scales above 100" or clumped on scales of a few arc minutes. 



Our models are compare d to the observations at 3.6 /im by IKashlinsky et al.l (l2007d. 

2005h a nd Cooray et al. ( 2007 ) in Figures US] and to observations at 1.6 /im from Thompson et al 
( l2007al ) in Figure [T71 For these observations, it is safe to treat them as "upper limits," as 
additional foreground contamination might still exist. At 3.6 /im, most of our predictions for 
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Fig. 16. — Our models for the angular p o wer sp ectra at 3.6 fim (halo+IGM) are compared 
with observations from iKashlinskv et al.l (12007a) (their Figure 1, lower panel, shown as the 



blue asterisks) and from ICooray et al.l (120071 ) (their Figure 1, images A, B, C, and D, with 



varying foreground galaxy cuts) shown as red diamonds. Most of our models lie beneath 
current observations. The mean intensity produced by Pop II stars with a Salpeter initial 
mass spectrum and /* = 0.5 is vl v = 15.1 nW m -2 sr _1 , which is over current observations. 
For Pop III stars with a Larson initial mass spectrum and /* = 0.01, vl v = 0.182 nW m -2 
sr -1 , which is allowed by observations (for more on the mean intensity, see section [9]). 
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Fig. 17. — O ur models for the angular power spectra (halo+IGM) are compared with obser- 
vations from iThompson et al.l (l2007al ) (for all sources deleted) at 1.6 /im, which are shown 
by the blue diamonds. Again, most of our models lie beneath current observations. As in the 
case at 3.6/im, the mean intensity from Pop II stars with a Salpeter initial mass spectrum 
and /* = 0.5 is high at vl v = 60.1 nW m" 2 sr^ 1 , while the mean intensity for Pop III stars 
with a Larson initial mass spectrum and /* = 0.01 is vl v = 0.802 nW m~ 2 sr" 1 . 
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Fig. 18. — Our models of the angular power spectrum (halos and the IGM) compared 
with the sens it ivitie s of upcoming CIBER missions (shown as the stepped blue lines) from 
Cooray et al.l ( 120091 ). CIBER will increase sensitivity of measured fluctuations, but still 
many of our models will lie beneath the detection limit. 



the angular power spectra are below the current observations, and are therefore still viable 
candidates. At 1.6 fim we see similar results. Therefore, it seems likely that early stars 
contribute at very low levels to the fluctuations in NIRB. Of course, other factors, such as 
the star formation time scale and the minimum redshift that star formation occurs at, z eil d, 
can also affect which models can agree with observations. 

Missions currently underway and future, more detailed experiments can make better ob- 
servations of the NIRB. AKARI (p reviously known as ASTRO-F) observed in 13 bands from 



2-160 jum (IMatsuhara et al.ll2008l ). The Cosmic Infrared Background Experiment (CIBER) 
will be able to obtain the power spectrum from 7" to 2 degrees. Combined with AKARI 
and Spitzer, fluctuations 100 times fainter than IRTS/DIRBE will be able to be observed. 
CIBER has two dual wide field imagers at 0.9 and 1.6 /im. An improved CIBER II will also 
measure fluctuations in four bands from 0.5 to 2.1 /im. This experiment will be pivotal to 



(Bock et al. 
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Coorav et al. 


2009) 



shown in F i gure [TK1 for both 0.9 and 1.6 /im (I and H-Band respectively) fjCooray et al.ll2009 



Bock et al 



20061 ). The sensitivity of CIBER will be much better than any of the current 



observations, but still many of our models lie beneath detection limits. 
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Table 3: Values of the star formation rate computed from our simulation, p\, in units of 
M Q yr _1 Mpc -3 . We have used the star formation timescale of tgp = 20 Myr. 



9. ADDITIONAL CONSTRAINTS FROM THE MEAN INTENSITY OF 
THE NEAR INFRARED BACKGROUND 

In addition to fluctuations, measurements have been taken of the mean intensity of the 
NIRB. Because these measurements rely on an accurate subtraction of the zodiacal light, 
measurements of the mean intensity of the NIRB are more difficult to perform. Currently, 
the interpretation of these measurements is still highly controversi al. Measurements of the 
excess in the NIRB (NIRBE) started out high (70 nW m" 2 sr' 1 ) jMatsumoto et all 12005 ) 
and have since declined. The most recent measurements are lower. iKashlinsky et al.l (l2007bl ) 
report that the mean intensity of the NIRBE must be greater than 1 nW m -2 sr" 1 to be 



consistent with fluctuations at 3.6 and 4.5 pm. iThompson et al. 



(I2007al) report a r esidual 



NIRBE of O.O^Qg at 1.1 and 1.6 pm. Fluctuations measured by ICooray et al.l ( 2007 ) imply 



that the mean NIRB cannot be much more than 0.5 nW m 2 sr 1 at 3.6 pm. Using these 
limits, can we put additional constrains on the first stars? 



We calculate the mean intensity of NIRB from ( lPeacocklll999l ) 

dzp([l + z]u, z) 



= 4vr J H(z)(l + z) 

where v is the observed frequency and p(u, z) is given by Eq. 
contained in p(u, z), is given by p*(z) = p*(z)/tsp(z), where 



P*( z ) 



(41) 

The star formation rate 
(42) 



where p\f°(z) is the mean mass density collapsed into halos taken from the simulation (which 
has the minimum halo mass of 2.2 x 10 9 M Q ), and is shown in the right panel of Figure |6j 
For the star formation timescale, we use tsF = 20 Myr, so that we can calculate the mean 
NIRB for models that are compatible with the WMAP data. The star formation rates with 
various star formation efficiencies are given in Table [3J 

We can now calculate the mean intensity of NIRB for our models with various stellar 
populations and values of /*. The value of / esc does not matter here because when calculating 
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Fig. 19. — Spectra of the NIRBE from various populations of stars over a redshift range of 
6 to 15. (Left panel) Changing initial mass spectrum and metallicity of the stars for a given 
star formation efficiency, /* = 0.01. (Right panel) Changing the star formation efficiency for 
Population III stars with a Salpeter initial mass spectrum. Models with high star formation 
efficiency, /* = 0.2, produce too high NIRB, and can be ruled out by the current upper limits 
from observations. Note that if we divided these curves by /*, they would become identical. 
In other words, these curves differ solely due to the varying values of /*. 



the mean intensity - it does not m atter where the photons are coming; from - the halo itself 
or the IGM surrounding the halo ([Fernandez fc Komatsull2006l ). The spectra of NIRB from 
various populations of stars (with varying mass, metallicity, and /*) over the redshift range 
of z = 6 — 15 (using simulation data up to the redshift 14.6) are shown in Figure [TU and 
their numerical values (integrated over 1 — 2 jttm) are tabulated in Table HI We also show 
the spectra of each radiation process in Figure [201 Finally, we show the mean intensity from 
two redshift bins, z = 6 — 10 and 10 — 15, in Figure [2~ll Lower redshift stars clearly dominate 
over stars at higher redshifts. This, combined with a sharp break due to the Lyman limit as 
well as a bump due to the Lyman-a line, may be used to constrain z cnd . 

Assuming that our equation for the star formation rate is accurate up to high redshifts, 
and using the parameters of this simulation, we can put constraints on the populations of first 
stars. If we take our u pper limit for the mean in tensity of the NIRBE to be 3 nW m -2 sr _1 
(the upper limit from iThompson et al.l (j2007al )). we can rule out most populations with 
high star formation efficiencies (/* = 0.2), unless star formation is constrained to only high 
redshifts. This is consistent with our fluctuation analysis - some of our models with high 
/* would produce angular power spectra above the levels observed. If the star formation 
efficiency is very low, say, /* = 0.001, then the mean background would be too small to 
detect. Most of the change in the amplitude of the NIRBE is from a change in the star 
formation efficiency while the metallicity and initial mass spectra of the stars affect the 
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Fig. 20. — Spectra of the NIRBE and how each component contributes to the overall inten- 
sity, over a redshift range of 6 to 15. 

shape of the spectra. Therefore, an accurate measurement of the mean NIRBE can give 
information on the star formation efficiency. Further constraints on the metallicity and mass 
may be possible with more precise observations in the future. 



10. PREDICTIONS FOR FRACTIONAL ANISOTROPY 

As we have seen, the magnitude of the predicted angular power spectrum depends on 
various parameters such as /*, t$F, /esc, and the initial mass spectrum. However, as the 
mean intensity also depends on these quantities, one may hope that the ratio of the power 
spectrum and the mean intensity squared would depend much less on these astrophysical 
parameters. 
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/. 


Redshift Range 


vl v 

Pop III Larson 


Pop III Salpeter 


Pop II Larson 


Pop II Salpeter 


0.2 


6-15 


32.2 


22.9 


43.2 


32.2 


0.01 




1.61 


1.15 


2.16 


1.61 


0.001 




0.161 


0.115 


0.216 


0.161 



Table 4: Values of the mean background intensity, vl v , in units of nW m~ 2 sr _1 for stars 
with different star formation efficiencies. The mean is calculated as an average of ul v over 1 
to 2 fim. 



Ignoring the IGM contribution and rewriting the halo contribution given by equa- 
tion fl37|) . we get 

c ( n b \ 2 r dz 

1 (4tt) 2 \ J *nJ J H{z)r\z){l + zY 

x \p\f\z) {[*(*) + (1 - / esc ) [Jff(z) + V\z) + P"( Z ) + }] 2 

X */(*=^) flta (*=^)- (43) 
By rewriting and averaging equation (141 p over a band, we get 

1 = 4^(%£) / H(z)(l + z) (44) 



x 



p h » l °{z) [F(z) + W(z) + V\z) + P"(z) + ^(z)} . (45) 



Therefore, in the ratio Ci/I 2 , /* and tgp (which is related to l a as / a oc I/^sf) cancel 
out exactly. The dependence on the initial mass spectrum, f(m), which determines l a via 
integral, nearly cancels out. However, the dependence on / esc does not cancel out: the power 
spectrum depends on / esc , whereas the mean intensity does not. Therefore, we conclude that 
the ratio depends primarily on / esc . In Figure 1221 we show the fractional anisotropy, 51 /I = 
+ l)Ci/ (27r/ 2 ), for various infrared bands. He re, I is the mean intensity a veraged over 



the bands defined in Table which are taken from lSterken fc Manfroidl (119921 ). We assume 
a rectangular bandpass. The upper curves are for / osc = 0.19, while the lower curves are for 
/esc = 1, which is consistent with the expectation: the ratio of the angular power spectrum 
to the mean intensity is lower for a higher / esc . We have checked that the ratio is nearly 
the same for different mass spectra for / osc = 0, in which case the dependence on l(z) nearly 
cancels out. 

Note that for / esc = the ratio, Ci/I 2 , may be regarded as an weighted average of 
b 2 eff (l/r)P lin (l/r). We find 51 /I = + 1)C , / /(2tt/ 2 ) w 10~ 2 with a weak dependence on 



-44- 




Fig. 21. — Spectra of the NIRBE for populations of stars over two redshift bins, z = 6 — 10 
and 10 - 15. 

I, i.e., 81 /I oc Z ' 25 . In other words, the expected fractional anisotropy of the near infrared 
background is of order a few percent for / esc = 0, and can be lower by a factor of a few for 



11. DISCUSSION AND CONCLUSIONS 

Any detection or non-detection of fluctuations in NIRB can give us information on stars 
forming at high redshifts, stars that could have helped to reionize the universe. The escape 
fraction of ionizing photons, the star formation efficiency, and the mass and metallicity of 
the stars can affect the amplitude and shape of the angular power spectrum of fluctuations 
in NIRB. 
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Band 


Center (microns) 


Waveband (microns) 


J 


1.25 


1.1-1.4 


H 


1.65 


1.5-1.8 


K 


2.2 


2.0-2.4 


L 


3.5 


3.0-4.0 


M 


4.8 


4.6-5.0 



Table 5: Band defi n itions used for infrared bands. These are given in Table 16.2 in 



Sterken fc Manfroidl (Il992h . 



We modeled the angular power spectrum from halos and the surrounding IGM by com- 
bining the analytic formulas for the luminosity of halos and the IGM with iV-body simula- 
tions coupled with radiative transfer for several different populations of stars. Shot noise is a 
major contributor to the angular power spectrum of halos at small scales, so it is important 
to include larger scales in observations to minimize the component of shot noise. 

The star formation efficiency has a significant effect on the amplitude of the angular 
power spectrum, with the amplitude of the angular power spectrum being proportional to 

For a fixed / T /tsF = fescf*Ni/ts? (i.e., for a given reionization history), a combination 
of parameters that maximize the star formation efficiency, /*, give the largest NIRB power 
spectrum; thus, the stars that are less massive and have more metals (smaller TVj), and are in 
halos with a lower escape fraction (smaller / esc ) will all increase the amplitude of the NIRB 
angular power spectrum. In general, the amplitude of the angular power spectrum of the 
halos (and the mean NIRBE) is mostly dominated by /*, while the angular power spectrum 
of the IGM can probe the ionization history though the factor f^/tgF- 

If we do not fix the reionization history, there are other parameters that can change 
the amplitude of the angular power spectrum significantly. The angular power spectrum is 
inversely proportional to the star formation time scale squared, C\ oc 1/tgF- This uncertainty 
in the star formation time scale can be directly related to our uncertainty in the mass to light 
ratio of galaxies, for igp oc M^jL^. As changes in tgp result in different reionization histories 
(for a given / 7 ), the other tracers of reionization, e.g., the electron-scattering optical depth 
measured by the WMAP satellite and the abundance of Lyman-a emitting galaxies, should 
help narrow down a range of magnitudes of the NIRB fluctuations that are consistent with 
what we already know about the cosmic reionization. 

The angular power spectrum of the IGM is typically a minor contributor to the overall 
fluctuations; however, the IGM contribution can be comparable to the halo contribution (the 
stellar contribution as well as the nebular contribution from within the halo), especially if 
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the escape fraction of ionizing photons from halos is high. In the limit that / csc is close to 
unity, we expect C/ GM /C; hal ° oc / e 2 sc , as Cf &l ° would be completely dominated by the stellar 
emission. One can even make the IGM contribution dominate over the halo contribution by 
increasing tgF, which suppresses the halo contribution as Cf al ° oc t$ F , and leads to a delay in 
the reionization. Yet, this would not change the IGM power spectrum significantly, as the 
IGM luminosity power spectrum saturates when the ionization fraction reaches X e ~ 0.5. 
Of course, we need to make sure that such a model can still complete the reionization by 
z ~ 6, and can reproduce the electron-scattering optical depth measured by WMAP. 

The redshift at which the formation of stars contributing to NIRB ends, z cn d, can also 
change the amplitude of the angular power spectrum significantly. Changing z en ^ affects not 
only the amplitude of the angular power spectrum, but also the shape. The attenuation of 
Lyman-a photons for the most part, does not affect the angular power spectrum of the halos 
greatly, but could affect the amplitude of the IGM by about a factor of 2. 



Previous estimates of the angular power spectrum of the NIRB by ICooray et al.l (120041 ) 
neglected to account for nonlinear bias. For our simulation with the minimum halo mass of 
2.2 x 10 9 M Q , the nonlinear bias is large enough to change the prediction for the shape of the 
angular power spect rum qualitatively: a turnover of l(l+l)Ci at / ~ 10 3 that was predicted by 



Cooray et al.l (120041 ) is not seen in our calculation, and the shape of the clustering component 
(i.e., minus the shot noise), is consistent with a pure power law, 1(1 + 1)C; oc I ' 5 . 

Note that our results for the shape of the angular power spectrum are valid for the 
minimum halo mass of M m - m = 2.2 x 10 9 M & . The non-linear bias would be smaller for 



smaller mass halos. For example, the simulation carried out by iTrac &: Cenl (120071 ) resolves 



halos down to a smaller mass, M min = 6 x 10 7 h^ 1 M Q , and would therefore find a smaller 
average bias. The halo bias is mass dependent, more massive halos being more strongly 
biased. The effective linear bias, 6 e ff,im; is the integral of the linear halo bias for a given 
mass, bi(Mh), times the mass function weighted by mass above a certain minimum mass, 
i.e., 6 e ff,iin = [J^ min dM h M h (dn h /dM h )bx{M h )]/[J^dM h M h {dn h /dM h )}. As there are 
many more halos at lower masses, lowering M min would result in a lower average bias. As 
the degree of the non-linear bias increases as the linear bias increases, one would find a 
smaller non-linear bias for a lower M min . Therefore, we might still see a turnover if these 
smaller halos are bright enough to contribute to the NIRB. 



However, the real situation would be more complex than the above picture. In llliev et al. 



(120071 ). we also performed reionization simulations which reso lve source ha l os do wn to this 



lower minimum mass of ~ 10 s M Q . There, however, unlike ITrac fc Cenl (120071 ). we took 
account of the fact that those small-mass (< 10 9 M ) halos are subject to Jeans filtering, 
and thus their star formation is suppressed if they reside within the ionized regions. The sup- 
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pression occurs disproportionately on the low-mass halos clustered around the high-density 
peaks (which are the first to be ionized). Therefore, the ultimate effect may not be as large 
as one might think by just including all halos down to ~ 10 8 M Q . While it is plausible that 
there may be a turn-over, it is also quite plausible that the location of the turn-over would 
be on a smaller scale than what would be predicted by the linear bias model. 

This c an, in principle, be studied using higher-resolution simulations like those in 



Iliev et al.l (120071 ) that include the Jeans filtering effect, but the simulation box size of 
35 h" 1 Mpc there is not quite large enough to give a reliable statistical measure of the 
large-scale structure and angular fluctuations in which we are interested. Towards that end, 
we have more recently performed a new set o f large-box, higher-re soluti on simulations tha t 



include the Jeans filtering effect, reported in IShapiro et al.l ( 120081 ) and llliev et al.l (j2008al ). 



We shall present results on the NIRB from these higher-resolution simulations elsewhere. In 
any case, the above consideration suggests that the shape of the angular power spectrum 
gives us important information about the nature of sources contributing to NIRB as well as 
the physics of cosmic reionization. 

Current observations seem to favor low levels of both the fluctuations and the mean 
NIRB due to the high- z (i.e., z > 7) sources. The current observations of the mean intensity 
of the NIRB seem to rule out high levels of /*, i.e., /* > 0.2. Most of our models for 
fluctuations still lie beneath the current observations of the fluctuations of the NIRB. The 
upcoming CIBER missions will improve the sensitivity of observations, but many of our 
models still lie below their sensitivity limits. Nevertheless, these new observations should be 
able to put tighter constraints on which high- 2 galaxy populations are allowed and which 
are ruled out. Given the lack of direct observational probes of high- z galaxy populations 
contributing to the cosmic reionization, the NIRB continues to offer invaluable information 
regarding the physics of cosmic reionization that is difficult to probe by other means. 
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A. Derivation of Angular Power Spectrum of NIRB Fluctuations 

The observed intensity (energy received per unit time, unit area, unit solid angle, and 
unit frequency) of the NIRB toward a direction on the sky n, I u (h), is related to the spa- 
tial distribution of the volume emissivity (luminosity em itted per frequency per comoving 
volume) at various redshifts, p(v, x, z), as (jPeacocklll999l p. 91) 



, . _ _c_ f p[u(l + z),hr(z),z] 
Mnj " AirJ H(z)(l + z) 



(Al) 



where r(z) = c dz'/H(z') is the comoving distance. The band-averaged intensity is then 
given by 



1(h) 



c f Qdvp[u(l + z),hr(z, 
dv i v \a) = — / dz- 
Air 



H(z)(l + z) 



dz 



(A2) 
(A3) 



4tt J "" H(z)(l + z) 2 

where v = u(l + z). Note that the denominator now contains (1 + z) 2 instead of (1 + z). 

Now, using the luminosity density integrated over bands, Pl(x, z ) = fvM+z) dv P^ v i x ' z )> 
we obtain 

— f dz Pl^W'*] 



^ - s j (A4) 

The spherical harmonic transform of 1(h), aim = J dh I(h)Y l ^ n (h), is then related to the 
three-dimensional Fourier transform of px(x, z), pl(\l,z) = j <i 3 x p^(x, ,z)e 4k ' x (the inverse 
transform is p(x, z) = J <i 3 k/(27r) 3 p L (\t, z)e~ lk ' yL ), as, 



a h 



dz 



4n(-i) 1 



d 3 k 



4tt J H(z)(1 + z) 2 
where we have used Rayleigh's formula: 

e -ik.ar W = A^^) l 3i\kr(z)]Y; m (k)Y lm (h). 



(A5) 



(A6) 



i in 



The angular power spectrum, C\ = (|a; m | 2 ), is then given by (() is the statistical ensem- 
ble average) 



C,= ^ 



dz 



Atx) J H(z)(1 + z) 2 J H(z')(1 + z') 2 



dz' 



- I k 2 dk P L (k,z)j l [kr(z)]j l [kr{z')] 

71 



(A7) 
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where we have used the definition of the luminosity-density power spectrum (Eq. (EZ 
the normalization of the spherical harmonics, J dkYJ m (k)Y^(k) = 1. 

Now, when / ^> 1, the integral within the square bracket can be approximated a: 



and 



k 2 dk Px(Mii[Mz)]iiM*')] 



5 [r(z) — r(z')] 



r 2 (z) 



Pl k 



I 



r(z) 



Therefore, 



C, 



where r 



-V 

AttJ 

c 

(4tt) 2 
r(z) and r' 



dz 



H(z)(l + z) 2 

dz 

H(z)r 2 (z)(l + z) 



dr' 



dz'/dr' 5(r — r') 



H(z')(l + z') 2 
I 



r 2 (z) 



Pi 



r(z) 



-Pl k 



r(z) 



(A9) 



r(z'), and we have used dz/dr = H(z)/c. This is Eq. (|3 



On the other hand, if one chooses to calculate C\ from a pair of /^(n) and /„'(n) instead 
of a pair of the band-averaged intensities i"(n), then (1 + z) A in the denominator of Eq. ( IA9I) 
becomes (1 + z) 2 : 



dz 



(4vr) 2 7 i/(z)r 2 (z)(l + z) 2 



P p [v{l + z),i/(l + z);k 



I 



r(z) 



(A10) 



where P p (u, v'\ fc , z) is the pow e r spec trum of the volume emissivity, p(u,z). Eq. (lAlOj) agrees 
with Eq. (10) of ICooray et al.l (120041 ) . Note that their j u is p(u)/(A7r) in our notation, which 
explains the absence of l/(47r) 2 i n their Eq. (10). However, this re sult does not agree with 



Eq. (3 ) of iKashlinsky et al.l (120041). which o r iginat es from Eq. (11) of lKashlinsky fc Odenwald 



( 120001 ). (See also Eq. (14) of IKashlinsky I ( 120051 ).) Kashlinsky et al.'s formula has (1 + z) 
in the denominator of Eq. (lAlOj) instead of (1 + z) 2 , and thus it misses one factor of 1/(1 + 
z). We were un a ble to trace the cause of this discrepancy. It is therefore possible that 
Kashlinsky et al.l ( 120041 ) over-estimated C\ by a factor of ~ 10. As they have not taken into 



account a short lifetime of massive stars in their calculations (i.e., they used Eq. (J7J) instead 
of Eq. (|6])), which yields another factor of ~ 100 in Ci, it is possible that they over-estimated 
Ci by a factor of ~ 10 3 . 



5 The exact integral of a product of two spherical Bessel functions is given by -| J k 2 dk ji(kr)ji(kr') = 
8(r — r')/r 2 . When a function, F{k), is a slowly-varying function of k compared to ji(kr)ji(kr'), which is a 
highly oscillating function for / ^> 1, we may obtain ^ J k 2 dk F(k)ji(kr)ji(kr') k, F(k — l/r)S(r — r')/r 2 . 
This is the so-called Limber's approximation. 
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Fig. 22. — Fractional anisotropy, 51 /I = \JT{1 + l)C//(27r/ 2 ), for different infrared bands, as 
labeled, versus the wavenumber I. This quantity depends primarily on the escape fraction, 
/es C . The upper curves are for / esc = 0.19, while the lower curves are for / esc = 1. Therefore, 
the expected fractional anisotropy of the near infrared background is of order of a few percent 
for / csc = 0, and can be lower by a factor of a few for / e sc = 1- Note that the dependence on 
/* and tsF cancels out exactly in 51 /I. 



